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Abstract. We present a new variational model for computing the electronic first- 
order density matrix of a crystalline material in presence of a local defect. A natural 
way to obtain variational discretizations of this model is to expand the difference Q 
between the density matrix of the defective crystal and the density matrix of the 
perfect crystal, in a basis of precomputed maximally localized Wannier functions of 
the reference perfect crystal. This approach can be used within any semi-empirical or 
Density Functional Theory framework. 
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Describing the electronic state of crystals with local defects is a major issue in solid- 
state physics, materials science and nano-electronics [Hill 13]. The first self-consistent 
electronic structure calculations for defective crystals were performed in the late 70', by 
means of nonlinear Green functions methods [H El [6]. In the 90', it became possible 
to solve the Kohn-Sham equations [7] for systems with several hundreds of electrons, 
and Green function methods were superseded by supercell methods [H |9]. However, 
supercell methods have several drawbacks. First, the defect interacts with its periodic 
images. Second, the supercell must have a neutral total charge, so that in the simulation 
of charged defects, an artificial charge distribution (a jellium for instance) needs to be 
introduced to counterbalance the charge of the defect. These two drawbacks may lead 
to large, uncontrolled errors in the estimation of the energy of the defect. In practice, 
ad hoc correction terms are introduced to account for these errors |10]. A refinement of 
the supercell approach, based on a more careful treatment of the Coulomb interaction, 
has also been proposed in [TT] . 

In a recent article [I2], we have used rigorous thermodynamic limit arguments to 
derive a variational model allowing to directly compute the modification of the electronic 
first order density matrix generated by a (neutral or charged) local defect, when the host 
crystal is an insulator (or a semi-conductor). This model has a structure similar to the 
Chaix-Iracane model in quantum electrodynamics [131 El- This similarity originates 
from formal analogies between the Fermi sea of a defective crystal and the Dirac sea 
in presence of atomic nuclei. For technical reasons, the reference model considered 
in [12] was the reduced Hartree-Fock model, or in other words, a Kohn-Sham model 
with fractional occupancies and exchange-correlation energy set to zero. 

The purpose of the present article is twofold. First, the extension of our model to a 
generic exchange-correlation functional is discussed. Second, a rigorous justification of 
the numerical method consisting in expanding the difference between the density matrix 
of the defective crystal and the density matrix of the perfect crystal, in a basis of well- 
chosen Wannier functions of the reference perfect crystal, is provided: this method can 
be seen as a variational approximation of our model. 

1. Derivation of the model 

We consider a generic Kohn-Sham model (or rather a generic extended Kohn-Sham 
model in which fractional occupancies are allowed) with exchange correlation energy 
functional E^^{p). For the sake of simplicity, we omit the spin variable. The ground 
state of a molecular system with nuclear charge density p™'^ and A/" electrons is obtained 
by solving 

inf {Efl^i^), < 7 < 1, Tr(7) = Af} , (1) 
Efl^^) = Tr (-1^^) - D{p--% p,) + \d{p,, p,) + E-(p,), (2) 
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where /^^(r) = 7(r, r) and where 

Jrs |r — r I 

is the Coulomb interaction. Still for simplicity, we detail the case of the Xa exchange- 
correlation functional 

Jm.3 

the extension to more accurate LDA functionals being straightforward. Likewise, 
replacing the all electron model considered here with a valence electron model with 
pseudopotentials does not bring any additional difficulty. 

The above model describes a finite system of Af electrons in the electrostatic field 
created by the density p™'^. Our goal is to describe an infinite crystalline material 
obtained in the thermodynamic limit JV —* oo. In fact we shall consider two such 
systems. The first one is the periodic crystal obtained when, in the thermodynamic 
limit, the nuclear density approaches the periodic nuclear distribution of the perfect 
crystal: 

p'^"^ - Pp:r^ (3) 

Pp"^ being a periodic distribution. The second system is the previous crystal in presence 
of a local defect: 

Typically, u describes nuclear vacancies, interstitial nuclei, or impurities together with 
possible local rearrangement of the nuclei of the host crystal in the vicinity of the defect. 
In the simple case of a reference perfect crystal with a single atom per unit cell 

„nuc \ ^ ^ r 
Pper = 2^ ^^R 

Ren 

where TZ is the Bravais lattice of the host crystal and 5r is the Dirac delta measure at 
R. If the defect consists in a impurity (the nucleus of charge z at R = being replaced 
with a nucleus of charge z'), the charge distribution u reads 

^ = z'Su(^o) - zSo+ z ((5r+u(r) - ^r) , 

Re7^\{o} 

where U is the displacement field of the nuclei generated by the relaxation of the crystal. 
It is therefore composed of nuclei of positive charges and of "ghost nuclei" of negative 
charges. In this article, we assume that u is given, and we focus on the calculation of 
the electronic density matrix. 

The form of the density matrix 7pgj. of the perfect crystal obtained in the 
thermodynamic limit ([3j) is well-known. The matrix 7°^^, is a solution to the self- 
consistent equation 

7per = X(-oo;.f](^Ser) (5) 
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Ker = -l^ + ^per-lc^aPll!\ (6) 

-A$pcr = 47r (p°,, - p^l^) , $pcr 7^-periodic. 

The notation P = X(-oo;ep](^) means that P is the spectral orthogonal projector of the 
self-adjoint operator A corresponding to filling all the energies up to the Fermi level ep 
(see for instance tl5j). In our case, ([5]) means that j^^^ is the spectral projector which 
fills all the energies of H^^^. up to the Fermi level ep, see Figure [H 

Zth band er {Z + l)st band 



7° 

' per 

Figure 1. Spectrum of H^^^. 

The density of the periodic Fermi sea is p^eri'^) = 7per(''5'')- Note that the system 
is locally neutral: 

/ r per / r per ' 

Jn Jn 

where is a reference unit cell, the Fermi level ep being chosen to ensure this equality. 
For the rest of the article, we assume that the host crystal is an insulator (or a semi- 
conductor), i.e. that there is a gap g = S+ — > between the highest occupied and 
the lowest virtual bands. Then the Fermi level can be any number S~ < ep < S"*". 

Now we consider the system obtained in the thermodynamic limit (jl]) when there 
is a defect u and derive a nonlinear variational model for it. We shall describe the 
variations of the Fermi sea with respect to the periodic state 7pgj,. The relevent variable 
therefore is 

Q ~ 1 ^ 7per 

where 7 is the density matrix of the defective Fermi sea. Notice that the constraint that 
7 is a density matrix (0 < 7 < 1) translates into — 7per < Q < 1 — 7per fo'^ the new 
variable Q. 

The energy of Q is by definition the difference of two infinite quantities: the energy 
of the state 7 and the energy of the periodic Fermi sea 7pgj.. Using ([2]), one obtains: 

S'^iQ) = Tt{HI,Q) - D{u,pq) + lDipQ,pQ) + 6-(pq) (7) 

where 

If we want to describe a defective crystal of electronic charge q {q electrons in excess 
with respect to the perfect crystal if g > 0, or —q holes if g < 0) interacting with the self- 
consistent Fermi sea in the presence of the defect, we have to consider the minimization 
principle 

E'^iq) = inf {S^iQ), -7°er < Q < 1 " 7Ser, MQ) = q} ■ (8) 
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We obtain in this way a model which apparently renders possible the direct calculation 
of the defective Fermi sea in presence of the nuclear charge defect u, when q electrons 
(or —q holes) are trapped by the defect. A globally neutral system would correspond to 
q = Jjg3 u but there is no obstacle in applying ([8]) to charged defects. 

Alternatively, one can, instead of imposing a priori the total charge q of the system 
(microcanonical viewpoint), rather fix the Fermi level ep G (S^,S+) (grand-canonical 
viewpoint). This amounts to considering the Legendre transform of ([8]): 

= inf {S-'iQ) - 6p Tr(Q), -7°,, < Q < 1 - 7p°cJ • (9) 
Any solution of ([H]) or ([9]) satisfies the SCF equation 

Q = X(-oo,.p) (^0) - 7Ser + 5, (10) 



where 

A 14 

and where < 5 < 1 is a finite-rank self-adjoint operator on L^(M^) such that 
Ran((5) C Ker(ifQ — ep). In the case of ([8]), the Fermi level ep is the Lagrange multiplier 
associated with the constraint Tr(Q) = q. The essential spectrum of Hq is the same as 
the one of H^^^ and is therefore composed of bands. On the other hand, the discrete 
spectrum of H^^^ is empty, while the discrete spectrum of Hq may contain isolated 
eigenvalues of finite multiplicities located below the essential spectrum and between the 
bands. Each filled (or unfilled) eigenvalue may correspond to electrons (or holes) which 
are trapped by the defect. 

Zth band ep (Z + l)st band 



7 = Q + 7per 
Figure 2. Spectrum of Hq. 



The SCF equation (ITOl) is equivalent to the usual Dyson equation, which is at the 
basis of Green function methods. 



2. Proper definition of the variational set 

The variational models ([S]) and may look similar to the usual Kohn-Sham models for 
molecules and perfect crystals. Their mathematical structure is however dramatically 
more complex. To design consistent numerical methods for solving and a deeper 
understanding of the mathematical setting is needed. 

The biggest issue with problems ([8]) and Qj is to properly define the variational set, 
that is the set of all Q's on which one has to minimize the energy functional £"{Q) or the 
free energy functional 8"{Q) — ep Tr(Q). For usual Kohn-Sham models, the variational 
set is very simple: it is the largest set of density matrices for which each term of the 
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energy functional is a well-defined number and the constraints are satisfied. This is the 
reason why it is not a problem to omit the precise definition of the variational set when 
dealing with usual Kohn-Sham models. For instance, the variational set for ([T]) is 

{7 I < 7 < 1, Tr(7) = Af, Tr(| V|7|V|) < 00} . (11) 

Let us recall (see [15] for instance) that if is a non-negative self-adjoint operator 
on L^(]R^) and if ('?/'«) jgN is an orthonormal basis of L^(M^), the series of non- negative 
numbers converges in M+ U {+C)o} towards a limit denoted by Tr(i?), 

which does not depend on the chosen basis. The operator B is said to be trace-class if 
Tt[B) < cxD. A bounded operator A on L^(M^) is trace-class if V A* A is trace-class. In 
this case, the scalar Tr(y4) = J2t^{'^i\^\^i) well-defined and does not depend on the 
chosen basis. On the other hand, if A is not trace-class, the series ^^^(V'il^l'^j) ^lay 
converge for one specific basis and diverge (or converge to a different limit) in another 
basis. 

The condition Tr(|V|7|V|) < 00 in ( fTTj) is a necessary and sufficient condition 
for each term of ([2]) being well-defined. In terms of Kohn-Sham orbitals, 
this conditions means that each orbital 0j is in the Sobolev space H^{M.^) = 
{0 G L2(R3) I V0 G (L2(M3))3}. 

The difficulty with the variational models ([8]) and ([9]) is that the variational set has 
not so simple a structure. It was shown in that an appropriate variational set is the 
convex set 

IC = {Q\ - 7°er <Q<1- ller, Tr(l + \V\)Q\1 + | V|) 

+ Tr(l + |V|)(Q++ - Q— )(1 + |V|) < 00}. 
In the above expression, we have used the notation 



with 

Q ~ 7perQ7per5 Q ~ IpcvQi^ ~ 7per)5 
= (1 - llMll^r. Q^^ = (1 - 7Ser)Q(l ' l^r) , 

corresponding to the decomposition 

L\M.^) = n-®n+, (12) 

where H- = 'j'^^^L'^ (M.^) and = (1 — ■y^^^) L"^ (M.^) are respectively the occupied and 
virtual spaces of the reference perfect crystal. 

Notice that when Q satisfies the constraint — 7per < Q < 1 — 7per5 o^i^ has Q~^~^ > 
and Q < 0. A remarkable point, proved in [12], is that the density pg of any operator 
Q G /C is a well-defined function which satisfies 

/ pI + DipQ^pq) < GO. 
Jr3 
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This shows that the electrostatic components of the energy £^^^(7) are well-defined and 
that so is the exchange-correlation contribution: as Pp^^. is periodic, continuous and 
positive on and as pq G L^(M^), the fifth term of which was not considered in 
[T2] is also well-defined. Finally, following [T6j, the generalized trace of an operator 
Q G /C is defined by 

Tr(Q) = Tr(Q++)+Tr(Q-), (13) 

and for any Q G /C, one sets 

Tr{H'^,,Q) = Tr([/7jJ++g++) +Tr([/7jJ-g-), 

where [H^er] [ifpgj.]++ are respectively the restrictions to the occupied and virtual 

spaces of the periodic Kohn-Sham hamiltonian of the perfect crystal. Note that H^^^. is 







[Kor]-- 1 








1 iKrV^ 



H = 

per 

The definition ffT3l) of the trace function is an extension of the standard trace function 



defined on the set of trace-class operators. Note that this extension depends of 7pgj, 
through the decomposition (|T2i) of the space. In the Quantum Electrodynamical 
model studied in [161 [13 HI IIS [H], minimizers are never trace-class (this property 
being related to renormalization). Whether or not the minimizers of (IHl) and iQ are 
trace-class still is an open question. 

To our knowledge, the variational interpretation of the ground state solutions of 
the self-consistent equation ( ITOl) as minimizers of the energy on the set /C with a 
constraint on the generalized trace (fT3l) . is new. This interpretation allows to rigorously 
justify the numerical method described in Section HI 



3. Interpretation in terms of Bogoliubov states 



The density matrix formalism used in the previous section can be reinterpreted in terms 
of Bogoliubov states, following [13] . 

Let 7 be an orthogonal projector acting on L^(]R^) such that Q = 7 — 7pcr ^ ^■ 
It can be proved [18] that there exists an orthonormal basis (0~)j>_Ar_ of 7i_ and an 
orthonormal basis {(pt)i>-N+ of TC+ such that in this basis 



Q 



\ 

















diag(-pi) 





diag(p-) 
















diag(p-) 





diag(pi) 



(14) 



with < Pi < I, Pi < 00, p[ = vPi(l — Pi). Notice that Q is a trace-class operator 
if and only if Ylt^ \fPi < 00 • Let us assume for simplicity that in equation (fTOl) . the 
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Fermi level ep is either empty or fully occupied. In this case, X(-oo,eF) (Hq) + 5 is an 
orthogonal projector, which implies that Q can be decomposed as in f|T4l) . It is important 
to mention that in this case, the generalized trace of Q is the integer — A^_. 

Formula (IT^ can be interpreted in terms of Bogoliubov states. The orbitals 
0^^^ , ■ ■ ■ , (p-i describe bound electrons in the virtual bands of the reference perfect 
crystal, while the orbitals (^I^ , ■ ■ ■ , 4>Zi represent bound holes in the occupied bands. 
Likewise, each pair {(f)f,(f)~) with i >0 and < < 1 is a virtual electron-hole pair, 
and (f)l and 0~ are the states of the corresponding Bogoliubov quasiparticles. The angle 
6i = asin(pj) is then called the Bogoliubov angle of the virtual pair. 

Formula fll4p can itself be rewritten in a second quantized form, using the Fock 
space built upon the decomposition f|T2l) . Let us introduce the A^-electron sector 
J-'^ := /\^ and the M-hole sector J^t^ := /\^^7i_. The electron-hole Fock space is 
defined as 

N,M>0 

We denote by a| the creation operator of an electron in the state (pf and by b] the 
creation operator of a hole in the state (p^ . In this formalism, the vacuum state 
f2o = 1 C?) 1 G (g) corresponds to the periodic Fermi sea of the perfect crystal, 
represented by the density matrix j^^^. in the usual Kohn-Sham description. We may 
also define the charge operator acting on the Fock space JF by 

i>-N+ i>-N- 

There is a special subclass of states in called Bogoliubov states [131 [IS 1201 [21]. 
Each Bogoliubov state G is completely characterized by its one-body density 
matrix 7, an orthogonal projector acting on L^(M^). Conversely, any projector 7 
gives rise to a Bogoliubov state under the Shale-Stinespring [221 [23] condition that 
Q = 1 ~ 7per is a Hilbert-Schmidt operator (which means Tr(Q^) < 00). The role 
of the Shale-Stinespring condition is to ensure that Q.y is a well-defined state in the 
same Fock space as the vacuum state Qq. Saying differently, this ensures that the Fock 
space representation associated with the splitting L^(M^) = 7L^(M^) © (1 — 7)L^(]R^) 
is equivalent to the one induced by ([I2D (i.e. ^^(^3) = 70^^l2(M3) © (1 _ 70^^1.2(^3)). 
Notice the Hilbert-Schmidt condition Tr(Q2) < cx) is satisfied for any Q = 7 — 7per in 
K,. Hence the variational set /C can be identified with a variational set of Bogoliubov 
states {fl^}^^jc in the Fock space J-'. 

The expression of the Bogoliubov state fl^ in the Fock space JF is given by [211 1221 [21] 




where Aj = tan(^i), and where c is a normalization constant. The above expression 
can be considered as the second-quantized formulation of (HM . It can then easily be 
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checked [12] that the charge of each Bogohubov state (counted relatively to that of 
the vacuum Qq) is actually given by (|T3|) : 

{n^\Q\n^) = Tr(Q++) + Tr(Q— ) = - N_ 

where Q = 7 - 7per- 



4. Variational approximation 



Let us now come to the discretization of problem ([8]). 

If one discretizes ([H]) in a local basis without taking care of the constraint Q G /C, 
there is a risk to obtain meaningless numerical results. On the other hand, selecting a 
basis set which respects the decomposition f|T2|) . will lead to a well-behaved variational 



approximation of ([H]) (the constraint Q G /C will be implicitly taken into acount). Let 



± 



be finite-dimensional subspaces of the occupied and virtual spaces H± of the reference 
perfect crystal. Consider the finite-dimensional subspace = V!^ © V]^ of L^(M^), the 
latter decomposition being the finite-dimensional counterpart of f|T2|) . Let (0i, ■ ■ ■ , 0m_) 
(resp. (0m_+i, • ■ • ,0AfJ) be an orthonormal basis of (resp. of V^). We denote for 
simplicity m+ := — m_. The approximation set for Q consists of the finite-rank 
operators 



(15) 



with eK,^ = [Q^ = [Q'^f, 0<T + Q'' <l}, where J is the x Nb block diagonal 
matrix 

lm,_ 





I 



The matrix of H^^j. in the basis (0j) is of the form 



H 



For Q of the form f|T5|) . it holds 







with 



and 



We then end up with the finite-dimensional optimization problem 
E,^(g) = inf {£'^m, e /C^ Tr(g'^) = q] 



(16) 
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which is a variational approximation of ([H]): 



As G with Tr(g'*) = g if and only if 

I + Q^ e {D = D'^ e R^^^ < D, Tr(D) = g + iV_} , 

problem ( |T6|) can be solved using relaxed constrained algorithms [25l [26] . 

The question is now to build spaces and that provide good approximations 
to ([8]) and ([9]). A natural choice is to use the maximally localized (generalized) Wannier 
functions [27] (MLWFs) of the reference perfect crystal. A very interesting feature of 
these basis functions is that they can be precalculated once and for all for a given host 
crystal, independently of the local defect under consideration. To construct Vf!:, one can 
select the maximally localized (generalized) Wannier functions of the occupied bands, 
that overlap with e.g. some ball Br^ of radius Rc centered on the nuclear charge defect. 
Note that due to the variational nature of the approximation scheme, enlarging the 
radius Rc systematically improves the quality of the approximation. To obtain a basis 
set for V"^ , one can select a number of active (unoccupied) bands using an energy cut-off 
and retain the maximally localized (generalized) Wannier functions of the active bands 
that overlap with the same ball B^^. The so-obtained basis set of the virtual space can 
be enriched by adding projected atomic orbitals of the atoms and ghost atoms involved 
in u (using the localized Wannier functions of the occupied bands to project out the TC- 
component of atomic orbitals preserves the locality of these orbitals). 

5. Numerical results 

In order to illustrate the efficiency of the variational approximation presented above, we 
take the example of a one-dimensional (ID) model with Yukawa interaction potential, 
for which the energy functional reads 



JR JR 

In the numerical examples reported below, the host crystal is Z-periodic and the nuclear 
density is a Dirac comb, i.e. 



with Z a positive integer. The values of the parameters {A = 10 and k = 5) have been 
chosen in such a way that the ground state kinetic and potential energies are of the 
same order of magnitude. 




with 
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Figure 3. Modulus of MLWFs associated witii tlie two occupied bands (left) and with 
the lowest two virtual bands (right). 



0^ 




Figure 4. Density pgh obtained with 28 MLWFs (line in red). The reference is a 
supercell calculation in a basis set of size 1224 (dashed line in blue). 



The nuclear local defect is taken of the form 



u= (Z - l)5o.25 - ZSq. 

This corresponds to moving one nucleus and lowering its charge by one unit. 

The first stage of the calculation consists in solving the cell problem. For simplicity, 
we use a uniform discretization of the Brillouin zone (— vr, tt], and a plane wave expansion 
of the crystalline orbitals. 

The second stage is the construction of MLWFs. For this purpose, we make use of 
an argument specific to the one-dimensional case [28]: the MLWFs associated with the 
spectral projector 7 are the eigenfunctions of the operator 7x7. One first constructs 
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A'^g mother MLWFs (taking 7 = 7pgj.)5 then Na mother MLWFs corresponding to the 
lowest Na virtual bands (taking for 7 the spectral projector associated with the lowest 
Na virtual bands). The so-obtained mother MLWFs are represented on Fig. [3l 

The third stage consists in constructing a basis set {4>j)i<j<Ni of N^, = Ny{Ne + Na) 
MLWFs by selecting the N^ translations of the {Ne + Na) mother MLWFs that are 
closest to the local defect, and in computing the first-order density matrix of the form 
f fTSj) which satisfies the constraints and minimizes the energy. The profile of the density 
Pqh obtained with Z = 2, Ne = 2, Na = 2 and Nb = 28 is displayed on Fig. HI It is 
compared with a reference supercell calculation with 1224 plane wave basis functions. 
A fairly good agreement is obtained with very few MLWFs. 

The implementation of our method in the Quantum Espresso suite of programs [29], 
in the true 3D Kohn-Sham setting, is work in progress [30]. 
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